Quantal corrections to mean-field dynamics including pairing 



Denis LacroisB and Danilo Gambacurta 
GANIL, CEA and IN2P3, Boite Postale 5027, 14076 Caen Cedex, France 



m 
o 

(N 



Sakir Ayik 

Physics Department, Tennessee Technological University, Cookeville, TN 38505, USA 

Extending the stochastic mean-field model by including pairing, an approach is proposed for 
describing evolutions of complex many-body systems in terms of an ensemble of Time-Dependent 
Hartree-Fock Bogoliubov trajectories which is determined by incorporating fluctuations in the initial 
state. Non-linear evolution of the initial fluctuations provides an approximate description of quantal 
correlations and fluctuations of collective observables. Since the initial-state fluctuations break the 
particle-number symmetry, the dynamical description in which pairing correlations play a crucial 
role is greatly improved as compare to the mean-field evolution. The approach is illustrated for a 
system of particles governed by a pairing Hamiltonian. 
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Under certain conditions, it is possible to provide an 
approximate description for quantal evolution of a sys- 
tem in terms of an ensemble of classical trajectories with 
proper choice of initial conditions. This aspect appears 
naturally in the path integral formulation of quantum dy- 
namics and has been recognized in refs. [HQ. I n recent 
years, this idea has been pushed forward to improved 
the mean-field description of many-body interacting sys- 
tems 0, 0). Mean- field description cannot describe es- 
sential quantal effects associated to collective motion and 
severely underestimates fluctuations of collective observ- 
ables. By considering an ensemble of mean- field trajec- 
tories with a specific choice of fluctuations (quantal zero- 
point and thermal) in the initial state, it is possible to 
overcome some of these shortcomings. Several applica- 
tions, especially in the nuclear physics context [5J, have 
shown that this approach can improve the mean-field de- 
scription by including important dissipative aspects in 
transport properties in heavy-ion collisions. More re- 
cently, we illustrated that such an approach, in addi- 
tion to fluctuations, can accurately tackle the problem of 
symmetry breakingclose to bifurcation point in collective 
energy landscape [6j. Up to now, the stochastic mean- 
field (SMF) approach with initial fluctuations has been 
developed starting from a Time-Dependent Hartrcc-Fock 
(TDHF) version of the mean-field. Nowadays, there are 
increasing interests in the treatment of pairing to describe 
evolution of strongly interacting Fermi liquids [Z - 11 em- 
ploying the Time-Dependent Hartree-Fock Bogoliubov 
(TDHFB) approach or its simplified BCS limit. While 
the TDHFB theory provides an important improvement 
beyond TDHF, it still suffers from the above quoted 
limitations: i.e. underestimation of quantum collective 
fluctuations and impossibility to spontaneously break 
symmetries. Due to successful description of the SMF 
approach in geometric symmetry breaking, it is rather 
tempting to introduce quantum fluctuations through ini- 



tial sampling to treat pairing where the U(l) symmetry 
breaking plays an important role. 

In this work, we present an extension of the SMF ap- 
proach by incorporating pairing into the description. In 
the standard TDHFB theory, evolution of the generalized 
density matrix lZ(t) is given by 



ih-ii(t) = [U{n),n{t)\ 



(i) 



lZ(t) contains both normal pij = (ajaj) and anomalous 

Kij = (cija,i) density matrix components. Here, {a\,ai) 
are creation and annihilation operators of a given single- 
particle basis. The quantity ri(7V) is the generalized 
mean-field Hamiltonian containing both the mean-field 
h[p] and the pairing field A[k] |l2j. Eq. (JTJ is gen- 
erally solved either starting from a density associated 
to a quasi-particle vacuum or a statistical ensemble of 
quasi-particles. It is then convenient to introduce asso- 
ciated quasi-particle creation operators, defined through 



7A 



J2i U* a ai + V* a a\ . In quasi-particle representation, 
normal density matrix is diagonal and expectation value 
of anomalous density matrix vanishes, 

Pa/3 = {Hill) = S a pfa, K a p = (7/37a) = 0, (2) 

where f a denotes the occupation numbers of quasi- 
particle states 12, [l3| and the expectation values arc 
taken with respect to the quasi-particle vacuum. As a 
result, lZ(t) takes diagonal form 

K(t) = \ W ^)) + l^»(*)> (1 - fa) (Va(t)\(3) 



1 In the following, Greek index will be reserved to quasi-particles 
while latin index will refers to particles. 
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where vectors \V a (t)) and | Wq, (£)) are the eigenvectors 
of lZ(t) whose expressions can be found in [12J. In fol- 
lowing, for notation compactness, we will use \S a (t)) 
for both types of states. Then, in the TDHFB dy- 
namics, we need to evolve state vectors according to 
ihd t \S a ) = H(lZ)\S a ) 7 while keeping the quasi-particle 
occupation fixed. In the stochastic extension of TD- 
HFB, quantal zero-point fluctuations (and thermal fluc- 
tuations) in the initial state are incorporated into the 
description in a similar manner to the SMF approach 
developed in [3[. Initial fluctuations are simulated by 
generating an ensemble of initial density matrices. Each 
member of the initial density matrix, labelled by (n) is 
written as, 



(4) 
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where summations over a and /3 run over a complete set 
of state vectors. The statistical properties of the elements 
of density matrix T^-^^o) is specified in terms of sta- 
tistical properties of the normal p^p(to) and anomalous 

K^g(io) density matrices. Elements of the normal and 
the anomalous density matrices are uncorrelated Gaus- 
sian random numbers with the mean values, 



P af} - o a pJ a , K a0 - U. 
and the second moments defined by, 
1 



(5) 



S & = §<W<W [/a(l - fp) + Ml - /a)] - ( 6 ) 

^S* = \s aa ,5 fW [U P + (1 - / a )(l - f p )] . (7) 



In these expressions, 5p( n > and Sk^ are the fluctuating 
parts of the density matrices, 

p( n ) = ~^n) + SpW , = ^ + (8) 

and notation X indicates the average over the initial en- 
semble. We note that in the SMF approach without pair- 
ing only fluctuations in the particle-hole channels are in- 
cluded. In the extension including pairing, additional 
fluctuations originating from particle-particle and hole- 
hole channels appears, Eq. ([7]). This is an important 
new aspect to explore initial conditions with non-zero 
anomalous densities. 

In the novel stochastic approach, an event of the en- 
semble of generalized density matrices is expressed as, 

n^{t) = £|si"^))<>o)<4 Tl) «|, (9) 

where states in an event evolve with its own self- 
consistent generalized Hamiltonian according to the TD- 
HFB equation of motion, ihd t \S ( a ] ) = H(ll in) )\si n) ). In 



the TDHFB framework, we can consider generalized op- 
erators that can create or annihilate two particles, 

Q = ^Ql/ata i + ^(Q2?a i a i + Q 3 2 fat a t).(10) 



The expectation value of such an operator in each event 
generated in the stochastic TDHFB approach reads 

Q (n \t) = Tr(QV n) W) + 29RTr(Q 2 V n) (t)). (11) 

Then, we determine the average evolution Q(t) by tak- 
ing the average over the ensemble generated in the sim- 
ulations. Since stochastic approach incorporates corre- 
lations beyond the mean-field description, the ensemble 
average value Q(t), in general, is different than the ex- 
pectation value of the operator in the standard TDHFB 
approach. Dispersion of the observable Q is calculated 
using the classical formula, 



4 



(t) = (QW(t)-Q(t)¥ 



(12) 



We note that, fluctuations (JBKZJ) are chosen to insure the 
initial ensemble average of the first and second moments 
equal to those obtained from the quasi-particle density 
matrix at zero or finite temperature. 

In order to illustrate the powerfulness of the approach, 
we consider a many-body system governed by a pairing 
Hamiltonian. In the model, there are K single-particle 
levels, labelled by i. Each level is associated with the en- 
ergy £j and has a degeneracy 2fij. It is assumed that each 
time a state i is present, its time-reversed states i is also 
present with same energy. We define number operators 
for each energy level as N t = Yip* ( a t,i a p,i + a p,i a P,i)- 
In addition, pair creation/annihilation operators Sf = 
-, \ t 

Sf ) arc introduced. The pairing 



Hamiltonian of the system is given by [14 



K 



K 



H = J2^i + J2 G ijSfSr 



(13) 



From the set of operators, one can construct three quasi- 
spin components Sf = (Sf+Sf)/2, Sf = (Sf - S^)/2i 
and Sf = (Ni — Oj)/2, which form a standard SU(2) al- 
gebra. The eigenstates of H can be determined by direct 
diagonalization (l5| giving access to the exact static and 
dynamical properties. 

Within TDHFB, evolutions of the expectation values 
(Sf(t), Sf(t), Sf(t)) of the quasi-spin components are de- 
termined by a set of coupled equations (h = 1), 

sm \ 

SUt) (14) 

sm ) 

where Af v (t) = ^2 k GikSf v (t). These equations of mo- 
tion are quite general and except that single-particle en- 
ergy are fixed in time, such equations can be found in 
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Fermi gas dynamics simulations on lattice with contact 
interaction [tJ or in nuclear physics context [lol 11|. 



In order to illustrate some limitations of TDHFB, let 
us suppose that the initial state of the system is de- 
scribed by a Slater determinant of single-particle wave 
functions. Since the Hamiltonian (|13[) contains two-body 
interactions, this state is not an eigenstate of the pairing 
Hamiltonian and it is expected to evolve in time. How- 
ever, it turns out that this state is a stationary solution 
of the TDHFB equation. Indeed, since the state is an 
eigenstate of the particle number, the expectation val- 
ues of quasi-spins vanish Sf(to) = Sf(t ) — 0. Hence, 
we deduce that all quasi-spin components remain con- 
stant in time. Therefore, the TDHFB approach is un- 
able to describe correlations that built up in time and 
leads to the departure from the independent particle ap- 
proach. This failure stems from the fact that a symmetry 
that exists initially, here the U(l) symmetry, cannot be 
spontaneously broken in the mean-field approach. In the 
stochastic TDHFB approach, the situation is different. 
Expectation values of quasi-spin components are still de- 
termined by the same set of equations (|14p. However, 
the initial condition is not determined in a determinis- 
tic manner, but is specified in terms of distributions of 
(S^ n \ S^ n \ S^ n ^) l= i.K- Following from the basic pos- 
tulate of the approach, © and ([7]) , the initial distribu- 
tions of quasi-spin components are uncorrelated random 
Gaussian numbers. Their mean values are determined 

-1/2 or -1/2 for 
Identify- 



by Sf n) = ST 1 ' = while St 
occupied and unoccupied levels respectively, 
ing the quasi-particle vacuum with a Slater determinant, 
variances of quasi-spin distributions are given by, 



?3/(") 



: z(n) 



gz(n) gz(n) 



: z(n) gy(n) 



Q x(n) Q y(n) 







and 
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are compared with stochastic simulations by generating 
iV ov t = 2 ■ 10 5 number of events. We obtain solutions of 
Eq. (fT4"|) by employing a Runge-Kutta 2 algorithm with 
a numerical time step At = 0.005/Ae. In Fig. [T] the 
quantity D(t), which illustrates the time-scale associated 
to the departure from the independent particle, 



D(t) = Y,ttini(t)(l - m(t)) 



(15) 



is plotted as a function of time for different coupling 
strength G/Ae. Here rij(t) are the single-particle occu- 
pation numbers. In this figure, symbols show the exact 
results and solid lines denote the simulations of stochastic 
TDHFB approach. In the stochastic TDHFB, the aver- 
age occupation numbers are determined by the following 

ensemble averages rii 




tune 



FIG. 1: (color online) Exact evolution of D(t) for G/Ae = 
0.05 (black triangles), 0.5 (blue squares) and 1 (red circles). 
The results obtained by averaging over TDHFB trajectories 
are shown respectively by black long dashed line, blue short 
dashed line and red solid line. Note that curves corresponding 
to G/Ae — 0.05 have been multiplied by 50. The time is 



presented in (Ae) unit. 



Consequently, for each level i, the z component of quasi- 
spin is initially a non-fluctuating quantity, while for each 
i, x and y components (S^ n \ Sf^) are specified by un- 
correlated real Gaussian distributions. Solution of Eq. 
(|14[) is determined by generating an ensemble of events 
according to the initial distributions of quasi-spin compo- 
nents. Solution is not a deterministic final state, but con- 
sists of a distribution of final states. Note that the mean 
number of particle is conserved event by event while the 
total energy is conserved in average. 

In the illustrations below, we consider K = 10 doubly- 
degenerated levels (ftj = 1) with constant level spacing 
Ae between adjacent levels, and constant pairing inter- 
action dk = G. A system of TV = 10 particles occupy 
initially the five lowest energy levels. The exact evolu- 
tion is obtained by decomposing the initial state in terms 
of eigenstates of the pairing Hamiltonian. The results 



The stochastic approach with pairing is able to incor- 
porate non-trivial dissipative effects due to the coupling 
between complex internal degrees of freedom. Incorpo- 
rating non-trivial initial fluctuations in particle-particle 
(pp) and hole-hole (hh) channels is able in all cases to 
describe short time dynamics very well. When the in- 
teraction strength increases, evolution of D{t) becomes 
more complex and tends to oscillate around an average 
value. The stochastic approach is accurate not only for 
short time behavior, but also gives a quite reasonable 
description of the asymptotic behavior although the evo- 
lutions appear more damped than the exact solutions. 
In the weak coupling regime G/Ae = 0.05, the exact dy- 
namics exhibits a periodic evolution with a rather short 
Poincare recurrence time that is perfectly reproduced by 
the average TDHFB evolution. We note that, with the 
present parameter value for the Hamiltonian, the BCS 
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threshold is at G cr /Ae ~ 0.3, i.e. above this coupling 
strength, a U(l) symmetry breaking mean- field solutions 
exists. It is important to recall that the standard TD- 
HFB gives D(t) = 0. Similarly, the previously proposed 
version of SMF approach that was not exploring the pos- 
sibility to break the particle number symmetry also leads 
to D(t) = 0. 

As another illustration, the two-body quantity a± (t) = 
^2ij{S^~S~)(t) is considered. Fig. [5] shows this quan- 
tity as a function of time for the exact and the stochas- 
tic TDHFB approach for different values of coupling 
strength G/Ae. Symbols indicate the exact evolu- 
tions and the solid lines are the results of the stochas- 
tic TDHFB simulations from the expression <r±(t) = 

Eij (sf n \t)Sf n) (t) + sf n \t)Sf n \t)) . Here also, we 



gross properties exact quantal dynamics. These illustra- 
tions indicate a strong support for the fact the stochas- 
tic extension of the mean-field without and with pairing 
provide a powerful tool for an approximate description 
of the correlations and quantal fluctuations of collective 
motion in normal and superfluid mesoscopic systems, be- 
yond the mean-field description. The TDHFB approach 
has recently became a widely employed tool to describe 
Fermi liquids appearing in nuclear, condensed matter or 
atomic physics. With rapid progress of computational 
powers, it is not anymore unreasonable to simulate thou- 
sands of independent TDHFB events, and consequently 
the stochastic approach proposed in this work may pro- 
vided as the next step to further progress in the descrip- 
tion of transport properties of such physical systems. 
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FIG. 2: (color online) Comparison of the exact and average 
evolution of u±(t) for different coupling strengths. Conven- 
tions are the same as in Fig. [T] Note that curves correspond- 
ing to G/Ae — 0.05 have been multiplied by a factor 5. 

observe the stochastic approach reproduces perfectly the 
exact evolution of quantal fluctuations during relatively 
short time scale for all coupling strengths, and provides 
a reasonable description of the gross properties (i.e. time 
average behavior) of the fluctuations over a long time 
scale for all coupling strengths 

In the present work, we propose a stochastic exten- 
sion of the mean-field dynamics including pairing corre- 
lations. In this approach, quantal zero-point (and possi- 
bly thermal) fluctuations in the initial state are incorpo- 
rated by simulating an ensemble of events by construct- 
ing a suitable distribution of the initial states. Observ- 
able quantities are determined by averaging over the gen- 
erated ensemble. The previous illustrations of the ap- 
proach without pairing [3, 6] and the present application 
including the pairing correlations show that the stochas- 
tic approach reproduces a reasonable description of the 
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